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ABSTRACT 

Ji. k — e model is proposed for wall bounded turbulent flows. In this model, the eddy 
viscosity is characterized by a turbulent velocity scale and a turbulent time scale. The 
thnp scale is bounded from below by the Kolmogorov time scale. The dissipation equation 
is reformulated using this time scale and no singularity exists at the wall. The damping 
function used in the eddy viscosity is chosen to be a function of Ry = instead of 

y+. Hence, the model could be used for flows with separation. The model constants used 
are the same as in the high Reynolds number standard k-e model. Thus, the proposed 
model will be also suitable for flows far from the wall. Turbulent channel flows at different 
Reynolds numbers and turbulent boundary layer flows with and without pressure gradient 
are calculated. Results show that the model predictions are in good agreement with direct 
numerical simulation and experimental data. 

I. INTRODUCTION 

In turbulence modeling, the Jfe-e model is the most widely used model in engineering 
calculations. The Standard Jfe - e Model 1,2 was devised for high Reynolds number turbulent 
flows and is traditionally used in conjunction with a wall function when it is applied to wall 
bounded turbulent flows. However universal wall functions do not exist in complex flows 
and it is thus necessary to develop a form of k — e model equations which can be integrated 
down to the wall. 

Jones and Launder 3 were the first to propose a low Reynolds number k-e model for near 
wall turbulence, which was then followed by a number of similar k-e models. A critical 
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evaluation of the pre-1985 models was made by Patel et al. 4 . More recently proposed 
models are found in Shih s and Lang and Shih 8 . Three major deficiencies can be pointed 
out about existing k - e models. (Some of the models may have only one or two of the 
three deficiencies.) First, a near wall pseudo-dissipation rate was introduced to remove the 
singularity in the dissipation equation at the wall. The definition of the near wall pseudo- 
dissipation rate was quite arbitrary. Second, the model constants were different from those 
of the Standard k — ( Model, making the near wall models less capable of handling flows 
containing both hi gh Reynolds number turbulence and near wall turbulence, which is often 
the case for a real flow situation. Patel et al. 4 put as the first criterion the ability of the 
near wall models to be able to predict turbulent free shear flows. Third, the variable y + is 
used in the damping function f M of the eddy viscosity formulae. Since the definition of y + 
involves ttr, the friction velocity, any model containing y + can not be used in flows with 
separation. 

In this paper, a new time scale based k - e model for near wall turbulence is proposed. 
In this model, k 1 ? 2 is chosen as the turbulent velocity scale. The time scale is bounded from 
below by the Kolmogorov time scale. When this time scale is used to reformulate the dissipa- 
tion equation, there is no singularity at the wall and the introduction of a pseudo- dissipation 
rate is avoided. The model constants are exactly the same as those in the Standard k - e 
Model, that ensures the performance of the model far from the wall. A damping function 
Is proposed as a function of Ry = ---y^ instead of y + . Thus, the present model can be used 
for separated flows. The performance of the model is then tested against turbulent channel 
flows at different Reynolds numbers and turbulent boundary layer flows at zero pressure 
gradient, favorable pressure gradient, adverse pressure gradient, and increasingly adverse 
pressure gradient. The results of the model prediction are compared with available data 
from direct numerical simulations and experiments. 

II. NEAR WALL * - e MODEL 

In turbulence modeling, the instantaneous quantities of an incompressible flow are de- 
composed into the mean and the fluctuating parts, i.e., u,- = U{ + Uj, p = P + p. The mean 
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field Ui satisfies the following continuity equation and Reynolds averaged Navier-Stokes 
equation: 

Ui,i = 0 ( 1 ) 

iii + UjUij = -- p Pi + vUiji- < «i«; >j ( 2 ) 

where the Reynolds stress term, — < t ii%ij >, must be modeled. 

In nn eddy viscosity model, one assumes mean field by 

— < XliUj >= Vr {Uij + Uj.i) ~ » ( 2 ) 

where v? is the eddy viscosity and Jb is the turbulent kinetic energy. From dimensional 
reasoning, the eddy viscosity is given by 

vj ~ tith, ( 4 ) 


where tit and l t are the turbulent velocity scale and turbulent length scale, respectively. In 
the framework of a high Reynolds number k — e model, u* ~ fc 1 / 2 , It ~ (which gives 
the turbulent time scale T t ~ *), and the resulting eddy viscosity is given by 

i 2 


vt = c. 


( 5 ) 


The transport equations for Jb and e are modeled as 

k + Ujkj = [(„ + J)y r < UiUj > Uij - € (6) 

c + Uj€j = [( v + - C\t < UiUj > U ijj - Clt-g 

= [(•' + — )«j]j + (~Cu < > Uij - C u *)/T t (7) 

Equations (5), (6), (7) with the following values of the constants = 0.09, C\ t = 1.44, 
C 2t = 1.92, < 7 * = 1.0, <7* = 1.3 are commonly referred to as the Standard k — ( Model 1,2 . 

A singularity would appear if the Standard k — c Model were applied down to the wall 
because of vanishing Jb at the wall, which renders the tune scale Tt (= k/c) in equation (7) 
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zero. However, as we shall show later, k/e can not represent the turbulent time scale near 
the wall. As the wall is approached, this commonly used time scale must be modified. 

The turbulent length scale is characterized by the size of the energy containing eddies. 
Near the wall, these eddies would have a size of 0(y). Following Hanjalic and Launder 7 , 
the turbulent velocity field has the following expansions near the wall: 

v! = <*iy + a 2 y 2 + — 
v' = b 2 y 2 + ... 

w' = ciy + c 2 y 2 + ... (8) 

where a*, b 2 , Ci are non-zero in general. Thus, as the wall is approached, both the turbulent 
length scale and the turbulent velocity scale approach zero. However, the turbulent time 
scale, which is given by the ratio of the length scale of the energy containing eddies to the 
turbulent velocity scale, approaches a non-zero value. We expect that this time scale must 
be the Kolmogorov time scale because viscous dissipation dominates near the wall. Thus, 
the turbulence time scale is given by k /e away from the wall and by the Kolmogorov time 
scale near the wall. Since Jb/e is much larger than the Kolmogorov time scale away from 
the wall and k/e vanishes near the wall due to the boundary condition on Jb, we can simply 
write 

T t = * + T fc (9) 

where 

T k = c k ty' 2 (10) 

is the Kolmogorov time scale and c k is a constant of order one. In the present investigation, 
c k = 1.0 is used. The constant c k was varied from 0.5 to 3.0 and the solutions were found 
to be quite insensitive to c k in this range. As more flow situations are tested, the value of 
c k could be optimized by fine tuning. 

Now the time scale given by equation (9) is bounded from below by the Kolmogorov 
time scale which is always positive. When this time scale is used in equation (7), there will 
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be no singularity at the wall. The argument on the time scale for near' wall turbulence and 
the method to remove the singularity in the dissipation equation at the wall are also used 
by other researchers, see Durbin 8 , for example. 

If we use Jfe 1 /* as the turbulent velocity scale and use the turbulent time scale as given 
by equation (9), the turbulent length scale is then the velocity scale times the time scale. 
Hence the eddy viscosity can be written as 

vr — fn^Tt ( 11 ) 


where /j, is the damping function which is used to account for the wall effect. 

The damping function in the present model is taken to be a function of Ry, which is 
defined by 


Ry = 


r > 


( 12 ) 


and takes the following form 

/^ = [1 - ezp(-aiRy - a 3 .R 3 - (13) 

where a a = 1.5xl0' 4 , a 3 = 5.0xl0~ 7 , o B = l.OxlO" 10 . These constants are devised by 
comparing the performance of the model prediction and the DNS data for turbulent channel 
flow at Re r — 180. 

Near the wall, the shear stress - < uv > should behave as 0(y 3 ), as derived from 
equation (8). Since Jfe is 0(y 2 ), we would require the damping function to have a near wall 
behavior of 0(y). From equation (12) and equation (13), it is seen that as y — ► 0, Ry — 1 ► 0 
as 0( y 1 ), which gives fy -* 0 as O(y). Thus, the near wall asymptotic behavior for the 
shear stress is satisfied. Far from the wall, Ry is large and / p — ♦ 1. The near wall eddy 
viscosity given by equation (11) then reduces to its counterpart for high Reynolds number 
flows given by equation (5). 

Since Ry instead of y + is used as the independent variable in the damping function, 
the model could be expected to be applicable in more general flow situations including 
flows with separation and reattachment. It should be pointed out that the use of y in 
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the damping function makes the model non-Galelian invariant. In addition, in complicated 
geometry situations such as comer flows, the meaning of y is ambiguous. However, in 
the authors’ view, these comments are more academic. In practical applications, the non- 
Galelian invariant will not be a problem through the use of a body fitted coordinate system 
with y denoting the direction normal to the wall. As for the comer flows, the ambiguity in 
y shows up only when it is very close to the comer, where k is very small, and the exact 
definition of y is not material. 

The other effect in near wall turbulence is the effect of the inhomogeneity of the mean 
field which introduces a secondary source term in the dissipation equation. Combining all 
the above, the dissipation equation is finally written as 

i 4- Uj€j = [(v + — + (— C\ t < UiUj > Uij — + E (14) 

<r* 

where the time scale T t is given by equation (9) and the secondary source term, E, is given 
by 

E = vvj Uij k U itjh . (15) 

This form was suggested by Jones and Launder 3 , and Shih 5 . The effect of E is confined to 
the buffer layer since away from the wall E becomes much smaller than the other terms in 
the dissipation equation. 


Equations (6), (11) (14) are the k - e equations proposed in this paper. The model 
constants, c M , C\ t i Cl <> are the same as those in the Stanford k - c Model. Away 

from the wall, the present model reduces to the Standard k — c Model. Thus, it is only 
necessary to assess the performance of the model for near wall turbulence. 


The boundary condition for e on the wall is determined by applying equation (6) (the 
k equation) at the wall, which gives 


f u ~ 


In this study, the following boundary condition for e, which is mathematically equivalent 
to the above but computationally much more robust, is used. 

(16) 


.dJb 1 / 2 ., 
= M-r-) 


dy 
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III. Numerical Aspects 


Bo undar y layer approximation is used in the calculations shown below. An semi-implicit 
finite difference scheme is used to solve the momentum equation and the transport equations 
for k and e. The coefficients for the convective terms are lagged one step in the marching 
direction and the source terms in the k and c equations are linearized in such a way that 
numerical stability is ensured. 

A variable grid spacing is used to resolve the sharp gradient near the wall. The grid 
distribution is controlled by Jy./tfy ;- 1 = a. Both a and the total number of the grid, N, 
are varied to ensure the grid independence of the numerical results. The marching step 
size, Sx, is also varied to ensure accuracy. In the test cases that follow, those parameters 
are changed such that the solution has a less than 1.0% error. Typically, the grids used 
are specified by N = 150 and a = 1.05. However, it is found that the solution is not very 
sensitive to the number of the grid points as long as there are two points in < 1. Fig. 1 
shows a calculation for the turbulent channel flow at Re r = 395 with N varying from 30 to 
150. It is seen that the results for different N are almost identical. 

IV. RESULTS AND DISCUSSIONS 

Turbulent channel flows at different Reynolds numbers and turbulent boundary layers 
with zero pressure gradient, favorable pressure gradient, adverse pressure gradient, and 
increasingly adverse pressure are calculated using the present model. The following shows 
the computational results compared with the available experimental data and data from 
direct numerical simulation. For some cases, the predictions of the Jones-Launder model 
and Chien’s model 9 are also shown. These two models are chosen because the Jones-Launder 
model is the first Jfc — e model for near wall turbulence while Chien’s model is known to 
perform quite well for turbulent boundary layer flows. 

Two dimensional fully developed channel flows were first chosen to validate the proposed 
model. These flows are attractive for model testing because they have self-similar solutions 
so that the initial conditions do not have to be accurately specified. These flows are very 
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simple and solutions can be found very efficiently; yet, the effects of the wall on turbulent 
shear flow are still present. In addition, DNS data providing detailed flow information is 
available for comparison. Computations are carried out for 2D fully developed turbulent 
channel flows at Re r = 180 and Re r = 395. Only the results for Re r — 395 are shown. 
Fig. 2 shows the profiles of the mean velocity, shear stress, turbulent kinetic energy and the 
dissipation. Both the dependent variable and the independent variables are represented in 
wall units by normalization through Ur and v. The predictions of the Jones-Launder model 
and Chien’s model are also shown. These predictions are compared with the DNS data. It 
is found that the present model predicts accurately the peak value of turbulent energy and 
the location of maximum value of the dissipation rate. (It should be at the wall.) Overall, 
the present model gives a better prediction. 

Like the 2D fully developed channel flows, turbulent boundary layer flows over a flat 
plate also give a self-similar solution. Thus, arbitrary profiles could be used as the initial 
conditions and the solution would develop into its similarity form. In our case, constant 
values are assigned to the velocity, turbulent energy and the dissipation rate. The exact 
values for the initial profiles are immaterial as long as the turbulent boundary gets generated. 
The boundary conditions used are that, at the free stream, the velocity reaches that of the 
freestream, while the turbulent energy and the dissipation rate are set to zero. Fig. 3 
shows the predicted velocity profile, turbulent energy, shear stress, and dissipation rate at 
Ret = 1410. The DNS data of Spalart 10 is also shown. Again, the predictions from the 
Jones-Launder model and Chien’s model are shown in the same figure for comparison. The 
results by the present model are found to be better. 

The computations were also carried out for larger Reynolds number cases and the com- 
parison is made between the model prediction and the experimental results of Wieghardt 
and W illmann 11 . Fig. 4 shows the results of the skin friction coefficient as a function of Reg. 
The skin friction by DNS at Reg = 1410 by Spalart 10 is also shown in this figure. Overall, 
the present model gives prediction of about 8% more than the experiment. It should be 
pointed out that for low Reynolds numbers situations, the skin friction from the experiment 
is lower than that from the DNS. For example, at Reg = 1410, C/(DNS) = 4.13 10 -3 , while 
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the ex per imen t gives C y (EXP ) = 3.97 10 3 , which is 4% lower than DNS. The present 
model predicts C/(Present) = 4.08 10" 3 . In Fig. 5, the velocity profile at Re 9 = 8900 is 
presented. It is seen that the present model prediction and the experiment agrees quite 
well. 

When the turbulent boundary layer is subject to a pressure gradient, the similarity 
solution ceases to exist. In this case, accurate description of the initial conditions for the 
velocity profile and the profiles of turbulence quantities (turbulent energy and dissipation 
rate) are very important. While experiment could provide the velocity profile at the up- 
stream location, information on the turbulence quantities is hardly available. In this study, 
the issue of the initial condition is dealt with in the following manner. We assume that 
the turbulent boundary layer develops under zero pressure gradient until it passes into 
the working section of the wind tunnel, where the experimental measurements are made. 
The connecting point between this virtual flat plate boundary layer and the real boundary 
layer with the pressure gradient is determined by the value of Re$ which is found by the 
experiment. 

The boundary conditions are specified in the same way as in the case of a flat plate 
boundary layer. At the wall, both the velocity and the turbulent energy are equal to zero 
while the dissipation rate is given by equation (16). At the free stream, zero values are 
assigned to the turbulent energy and its dissipation rate. The mean velocity approaches 
that of the free stream, which is determined by the experiment and is related to the pressure 
gradient of the flow. 

In the present study, the turbulent boundary layer with favorable pressure gradient 
studied by Herring and Norbury 11 , the turbulent boundary layer under adverse pressure 
gradient studied by Bradshaw 11 , and the turbulent boundary layer under increasingly ad- 
verse pressure gradient studied by Samuel and Joubert 13 are calculated, respectively. The 
first two are the test cases for the 1968 Stamford Conference on Turbulent Boundary Layers 11 
and the last one is a test case in the 1981-82 Stanford Conference on Complex Turbulent 
Flows 13 . The results of the calculation are shown below and the comparison is made with 
the experimental data in each of the cases. 
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Herring and Norbury 11 studied the development of a turbulent boundary layer under 
a favorable pressure gradient. At the first point of the working section of the experiment, 
Re 0 = 3400. Thus, a flat plate boundary at Re$ = 3400 is used to provide the initial 
conditions. With the initial conditions given, the calculation of the boundary layer is then 
carried out downstream. The result for the skin friction coefficient is shown in Fig. 6, and 
the result for the mean velocity at x = 4ft is shown in Fig. 7. The distances, x and y, are 
in physical units while the skin friction and the velocity are normalized by the freestream 
velocity at the streamwise location under consideration. It is seen that both the velocity 
and the skin friction are well predicted. 

The Bradshaw 11 experiment simulated in this paper is a flow where the turbulent bound- 
ary layer develops under an adverse pressure gradient. The first point where the velocity 
profile (and hence Ref) is known is at x = 2ft, where Ret = 7810. The experiment clearly 
indicates that the pressure gradient is quite strong at x =2ft. (The strongest over all the 
working section for this experiment.) Thus, the specification of the initial conditions based 
on a fiat plate boundary layer may not be very accurate. This is evidenced by Fig. 8a which 
shows the velocity profile from the experiment and the velocity profile from the flat plate 
boundary calculation. The difference between the two is quite substantial. Our prediction 
then carries some error with it due to the specification of the initial conditions; however, 
we hoped that this error would decrease as the calculation marches downstream. This is 
indeed found to be true. Fig. 8b shows the velocity profiles at x= 4ft. It is clear that 
the model gives quite good results even though the initial conditions are not accurately 
prescribed. The distribution of the skin friction coefficient is shown in Fig. 9. As in the 
case of the favorable pressure gradient, both x and y are shown in physical units, while the 
skin friction and the velocity are normalized by the freestream velocity at the streamwise 
location under consideration. 

In the Samuel and Joubert flow, the turbulent boundary layer develops under increas- 
ingly adverse pressure gradient. During the 1981-82 Stanford Conference on Complex Tur- 
bulent Flows 13 , it was found that this flow is very difficult to predict. Since then, it has 
become a strong test case for turbulence models. 
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The initial conditions are specified at x = 0.855m, where Re t = 5470. First, the pre- 
dicted skin friction is shown in Fig. 10 with the experimental values. The variation of the 
skin friction coefficient with x is well-captured. However, the predicted value is higher than 
the experiment result by a constant amount of 3.0 X 10 which amounts to about 11% 
at * = 0.855. It should be mentioned that at * = 0.855 where the pressure gradient is very 
small, one would expect that the experimental result for the skin friction coefficient would 
be the same as that for the flat plate boundary layer at the same value of Re t , found by the 
experiment of Wieghardt and Willmann, for example. However, it is seen that the result for 
this case is smaller than that of Wieghardt and Willmann by 1.4 X 10 -4 , which amounts 
to 5%. The velocity profile at x = 1.76m is shown in Fig. 11. The turbulence quantities 
(shear stress and the turbulent energy) at x = 1.79m are shown in Fig. 12 and Fig. 13, 
respectively. 


V. CONCLUSION 

We have presented a J fc — c model for wall bounded flows which is free from the three 
deficiencies we pointed out in the introduction. First, the proposed model uses the same 
set of model constants as that used in the Standard k — t Model and away from the wall the 
proposed model will reduce to the Standard k - e Model. Thus, the proposed model would 
be applicable in both near wall turbulence and high Reynolds number turbulence. Second, 
the proposed model uses a time scale which has the Kolmogorov time scale as its lower 
bound. By using this time scale to reformulate the dissipation equation, the singularity in 
the standard dissipation equation is removed as the wall is approached and the equation can 
be integrated to the wall. This renders the introduction of pseudo-dissipation unnecessary. 
Third, the proposed model uses Ry instead of y + as its independent variable in the damping 
function. This allows the model to be used in more complicated flow situations, flows with 
separation, for example. 

Turbulent channel flows at different Reynolds numbers and turbulent boundary layers 
with/ without pressure gradient are calculated using the present model. At low Reynolds 
number, the comparison between the DNS data and the present model is found to be excel- 
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lent. At higher Reynolds numbers, the velocity profiles are well-predicted. The predicted 
skin friction coefficient is within 8% of that of the experiment for a flat plate boundary 
layer, and is in an excellent agreement with the boundary layer under favorable pressure 
gradient (Herring and Norbury flow). For adverse pressure gradient cases, the predicted 
skin friction coefficient is less than 18% higher for the Bradshaw flow and for the Samuel 
and Joubert flow. 

It should be mentioned that the model is computationally robust. Arbitrary initial pro- 
files can be used for turbulent channel flows and flat plate boundary layers when similarity 
solutions exist. The predicted solution is also found to be quite insensitive to the n umb er 
of grid points. 
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Fig. 1: Mean velocity profile found by different number of grid points. 
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Fig. 2: Turbulent channel flow at Re r — 395 
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Fig. 3: Turbulent flat plate boundary layer at Re$ = 1410 
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Fig. 4: Skin friction distribution for turbulent flat plate boundary layer. 



Fig. 5: Velocity profile for turbulent flat plate boundary layer at Re t = 8900 
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Fig. 6: Skin friction distribution for Herring and Norbury flow. 
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Fig. 7: Velocity profile for Herring and Norbury flow at x = 4ft. 
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Fig. 9: Skin friction distribution for Bradshaw flow. 
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Fig. 10: Skin friction distribution for Samuel and Joubert flow. 
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Fig. 11: Velocity profile for Samuel and Joubert flow at x = 1.76m. 
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Fig. 12: Turbulent energy for Samuel and Joubert flow at x = 1.79m. 
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Fig. 13: Turbulent shear stress for Samuel and Joubert flow at i = 1.79m. 
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